\input{euler.tex}

\begin{document}

\problem[360]{Scary Sphere}

Given two points $(x_1, y_1, z_1)$ and $(x_2, y_2, z_2)$ in three dimensional space, the \emph{Manhattan distance} between those points is defined as $|x_1-x_2|+|y_1-y_2|+|z_1-z_2|$. 

Let $C(r)$ be a sphere with radius $r$ and center in the origin $O(0, 0, 0)$. Let $I(r)$ be the set of all points with integer coordinates on the surface of $C(r)$. Let $S(r)$ be the sum of the Manhattan distances of all elements of $I(r)$ to the origin $O$.
 
For example, $S(45) = 34518$. 

Find $S(10^{10})$.

\solution

There are three main methods to deal with this problem. The first method is a generic algorithm that generates lattice points on the sphere using Pythagorean quadruple formulas. It can be applied to any sphere lattice related problem. The second method makes use of the symmetry of Manhattan distance and sphere, and computes the result from the sum of squares function. This method is restricted to computing the Manhattan distance, but allows spheres of any integral radius. The third method is specifically tailored to deal with spheres of radius $2^n 5^n$, and presents a recursive formula to generate the lattice points quickly. This method provides some insight into the layout of the lattice points.

At present, only the first method is implemented, and it is described below. A brute-force method is also implemented as a benchmark for small scale problems.

Common to all methods, it's useful to note the fact that in any primitive Pythagorean quadruple, $a^2+b^2+c^2=d^2$ where $\gcd(a,b,c,d)=1$, $d$ must be odd. This means that for any lattice point $(x,y,z)$ on a sphere of radius $2^n r'$, $x, y, z$ must be divisible by $2^n$. Hence $S(2^n r') = 2^n S(r')$.

\method{I}

Without loss of generality, we can just consider the integer points on the surface of the sphere in the first octant, i.e. $x, y, z \ge 0$. There are three types of integer points:

I. Points on the pole, e.g. $(0,0,z)$. Each point of this type in the first octant generates six points on the sphere (including itself).

II. Points on the main orbit, e.g. $(0,y,z)$. Each point of this type in the first octant generates twelve points on the sphere (including itself).

III. Points within the first octant and not on any axes. Each point of this type in the first octant generates eight points on the sphere (including itself).

To enumerate the lattice points $(x, y, z)$ on the sphere of radius $r$, note that $(x, y, z, r)$ form a Pythagorean quadruple, which can be generated by $(ka, kb, kc, kd)$ where
\begin{align}
a &= m^2 + n^2 - p^2 - q^2 \notag \\
b &= 2(mq + np) \notag \\
c &= 2(nq - mp) \notag \\
d &= m^2 + n^2 + p^2 + q^2 \notag
\end{align}
where $\gcd(m,n,p,q)=1$ and $m+n+p+q \equiv 1$ (mod 2). Hence, for each odd divisor $d$ of $r$, we can enumerate all $(m,n,p,q)$ quadruples that satisfy the above conditions to find out the lattice points.

To enumerate all $(m, n, p, q)$ quadruples such that $m^2 + n^2 + p^2 + q^2 = d$, we require that $m \le n \le p \le q$ for the enumeration and then permute the result. To reduce running time, we first compute and store all possible sums of squares of the two smaller terms, $m^2 + n^2$. Note that
\begin{align}
d \ge 4 m^2 &\Rightarrow m \le \sqrt{d} / 2 \notag \\
d \ge m^2 + 3 n^2 &\Rightarrow n \le \sqrt{(d-m^2)/3} , \notag
\end{align}
so the total number of terms stored is
\begin{align}
&\sum_{m=0}^{\lfloor \sqrt{d} / 2 \rfloor} \sum_{n=m}^{\lfloor \sqrt{(d-m^2)/3} \rfloor} 1 \notag \\
&= \sum_{m=0}^{\lfloor \sqrt{d} / 2 \rfloor} \left( \left\lfloor \sqrt{\frac{d-m^2}{3}} \right\rfloor -m+1 \right) \notag \\
&\le \sqrt{\frac{d}{3}}+1 + \int_0^{\sqrt{d}/2} \left(\sqrt{\frac{d-x^2}{3}}-x+1\right) \, dx \notag \\
&= \sqrt{\frac{d}{3}}+1 + \frac{1}{\sqrt{3}} \left. \left( \frac12 x \sqrt{d-x^2} + \frac{d}{2} \arctan \frac{x}{\sqrt{d-x^2}} \right) \right|_{0}^{\sqrt{d}/2} - \frac{d}{8} + \frac{\sqrt{d}}{2} \notag \\
&= \frac{\pi}{12 \sqrt{3}} d + \left(\frac{1}{\sqrt{3}} + \frac12 \right) \sqrt{d} + 1 \notag \\
&\le 0.15115 d + 1.07736 \sqrt{d} + 1 . \notag
\end{align}

\complexity

Let $d$ be the largest odd divisor of $r$.

Time complexity: $\BigO(d \ln d)$

Space complexity: $\BigO(d)$

\method{II}

Not implemented.

%Since the Manhattan distance to the origin is just the sum of the absolute value of each coordinate, the total distance is the sum of the absolute coordinate value of all lattice points. Suppose $z > 0$ and $(x,y,z)$ is a lattice point, then all six points below are lattice points
%\[
%(x,y,\pm z), (\pm z,x,y), (y, \pm z, x).
%\]
%To ensure that these points are distinct, we require $z$ to be odd. This can be shown as below. Since $r$ is odd by construction, we have $r^2 \equiv 1$ (mod 4). Since $x^2+y^2+z^2 = r^2$, it follows that exactly one of $x, y, z$ is odd. (Otherwise, if all are odd, then the right hand side would be equal to 3 mod 4.) So without loss of generality, we can assume $z$ is odd.
%
%By symmetry, each $z$ is added six times. So the total distance is $6z$ times the number of solutions for each $z$.

% Note that
% \[
% x^2 + y^2 = r^2 - z^2 = (r+z)(r-z)
% \]
% For a given $z$, the number of solutions to the above equation can be denoted as $r_2(r^2 - z^2)$, where $r_2(\cdot)$ is the \emph{sum of squares function}. It can be computed by factorizing $(r+z)(r-z)$ as (if possible)
% \[
% (r+z)(r-z) = 2^{a_0} p_1^{2a_1} \cdots p_r^{2a_r} q_1^{b_1} \cdots q_s^{b_s}
% \]
% where the $p_i$s are primes of the form $4k+3$ and the $q_i$s are primes of the form $4k+1$. If such factorization is possible, then 
% \[
% r_2(\cdot) = 4 (b_1+1) \cdots (b_s + 1) .
% \]

\method{III}

Not implemented.

\answer

878825614395267072

\reference

http://en.wikipedia.org/wiki/Pythagorean\_triple

http://en.wikipedia.org/wiki/Pythagorean\_quadruple

http://integral-table.com

%http://mathworld.wolfram.com/SumofSquaresFunction.html

%http://www.fq.math.ca/Scanned/34-2/oliverio.pdf

%http://ebooks.library.cornell.edu/cgi/t/text/pageviewer-idx?c=math\&cc=math\&idno=05150001\&node=05150001\%3A4\&view=image\&seq=45

%http://mathworld.wolfram.com/CircleLatticePoints.html

\end{document} 